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I.  EXECUTIVE  SUMMARY 

The  main  objectives  of  this  project  were:  to  advance  the  understanding  of  the  generation, 
propagation  and  breakdown  of  atmospheric  gravity  waves  owing  to  mountains;  to  develop 
improved  theoretical  models  for  forecasting  stratospheric  gravity  waves  under  realistic  flow 
conditions;  and  to  examine  the  role  of  gravity  waves  in  the  production  of  turbulence  under 
various  flow  conditions. 

Towards  reaching  these  goals,  the  following  specific  problems  have  been  studied: 

(i)  As  indicated  by  balloon  measurements,  atmospheric  buoyancy  profiles,  apart  from  a 
sharp  increase  (roughly  by  a  factor  of  2)  at  the  tropopause,  often  feature  appreciable  oscilla¬ 
tions  (typical  wave  length  1-2  km)  with  altitude.  Based  on  numerical  solutions  of  the  Euler 
equations,  we  find  that  such  short-scale  oscillatory  variations  of  the  background  buoyancy 
frequency,  which  usually  are  ignored  in  theoretical  models,  can  have  a  profound  effect  on 
the  generation  of  mountain  waves  owing  to  a  resonance  mechanism  that  comes  into  play  at 
certain  wind  speeds  depending  on  the  dominant  oscillation  wave  length.  As  suggested  by 
our  earlier  asymptotic  results  for  small-amplitude  sinusoidal  oscillations,  comprehensive  nu¬ 
merical  studies  for  more  realistic  flow  conditions  confirm  that,  under  resonant  conditions, 
the  induced  gravity-wave  activity  is  significantly  increased,  often  resulting  in  wave  breaking 
above  and  upstream  of  the  mountain,  similarly  to  resonant  flow  of  finite  depth  over  topog¬ 
raphy.  This  piece  of  work  forms  part  of  the  doctoral  thesis  of  I.  Skopovi  (completed,  May 
2006).  For  technical  details,  see  Sec.  A  below. 

(ii)  The  effects  of  unsteadiness  in  the  wind  have  received  relatively  little  attention  in  prior 
studies  of  mountain  gravity  waves.  To  fill  this  gap,  the  effect  of  temporal  variations  in  the 
wind  velocity  was  studied.  Motivated  by  field  observations,  which  typically  exhibit  strong 
oscillatory  fluctuations  about  a  steady  mean,  we  considered  the  case  of  a  time-harmonic  com¬ 
ponent  superposed  on  a  constant  wind.  We  find  that,  for  relatively  short  (on  the  order  of  min.) 
wind  fluctuations,  the  unsteady  wind  component  can  be  safely  averaged  out.  On  the  other 
hand,  for  larger-period  (hours  or  days)  fluctuations,  the  effect  of  unsteadiness  is  very  signifi¬ 
cant.  In  particular,  steady  flow  states,  predicted  on  the  assumption  of  uniform  wind,  may  not 
be  attainable  owing  to  transient  wave  breaking.  This  forms  the  second  part  of  I.  Skopovi’s 
doctoral  thesis  (May  2006).  For  technical  details,  see  Sec.  B  below. 

(iii)  Recent  balloon  measurements  in  the  course  of  the  French  gravity  wave  campaign  in¬ 
dicate  that  low-frequency  (hydrostatic)  gravity  waves  are  ubiquitous  in  the  atmosphere.  These 
are  large-scale  disturbances  (horizontal  length  of  hundreds  of  km)  for  which  the  effects  of  the 
Earth’s  rotation  cannot  be  neglected.  In  an  effort  to  understand  how  energy  gets  transferred  to 
such  low-frequency  motions,  we  studied  the  nonlinear  evolution  of  a  non-hydrostatic  wave- 
train  as  it  propagates  upwards  in  the  atmosphere.  We  find  that  the  envelope  of  the  wavetrain 
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can  interact  resonantly  with  low-frequency  gravity-inertial  wave  modes.  This  resonant  inter¬ 
action  causes  radiation  of  energy  away  from  the  gravity  wavetrain,  forming  a  trailing  wake 
of  gravity-inertial  waves.  Moreover,  the  gravity  wavetrain  may  experience  rapid  wavenum¬ 
ber  variations  in  certain  locations,  where  eventually  caustics  form,  signaling  the  onset  of  a 
breakdown  process.  This  is  joint  work  with  Dr.  Ali  Tabaei,  who  was  supported  by  the  grant 
as  a  post-doctoral  research  associate.  Techinical  details  can  be  found  in  a  forthcoming  article 
in  Studies  in  Applied  Mathematics  (Tabaei  &  Akylas  2007). 

-  Personnel  supported 

Professor  T.R.  Akylas,  Dr.  A.  Tabaei  (post-doctoral  Research  Associate),  I.  Skopovi  (grad¬ 
uate  Research  Assistant) 

-  Publications 

Skopovi,  I.  &  Akylas,  T.R.,  “The  Role  of  Buoyancy-Frequency  Oscillations  in  the  Gen¬ 
eration  of  Mountain  Gravity  Waves”,  Theoretical  and  Computational  Fluid  Dynamics,  in 
press. 

Skopovi,  I.  &  Akylas,  T.R.,  “Effects  of  Temporal  Variations  in  Wind  Velocity  on  the  Gen¬ 
eration  of  Mountain  Gravity  Waves”,  Theoretical  and  Computational  Fluid  Dynamics,  to 
be  submitted. 

Tabaei,  A.  &  Akylas,  T.R.,  “Resonant  Long-Short  Wave  Interactions  in  an  Unbounded 
Rotating  Stratified  Fluid”,  Studies  in  Applied  Mathematics,  in  press. 

Akylas,  T.R.,  “Solitary  Waves  in  Rotating  Fluids”,  Chapter  5  in  Solitary  Waves  in  Fluids 
(ed.  R.H.J.  Grimshaw)  WIT  Press  2007. 

Akylas,  T.R.,  Grimshaw,  R.H.J.,  Clarke,  S.R.  &  Tabaei,  A.,  “Reflecting  Tidal  Wave  Beams 
and  Local  Generation  of  Solitary  Waves  in  the  Ocean  Thermocline”,  Journal  of  Fluid 
Mechanics,  under  review. 

Skopovi,  I.,  “The  Role  of  Background  Flow  Variations  in  Stratified  Flows  Over  Topogra¬ 
phy”,  PhD  Thesis,  May  2006. 


1 


m 


3 


II.  TECHNICAL  DESCRIPTION 

A.  Effects  of  Buoyancy-Frequency  Oscillation  on  the  Generation  of  Mountain  Waves 

The  flow  of  wind  over  mountain  ranges  is  believed  to  be  one  of  the  primary  generation 
mechanisms  of  gravity  waves  in  the  atmosphere.  Typically,  mountain  waves  develop  instabil¬ 
ities  as  they  propagate  upwards  in  the  atmosphere  and  eventually  breakdown  into  turbulence. 
There  is  now  increasing  evidence  (Eckermann  and  Preusse  [11])  that  breaking  mountain 
waves  in  fact  contribute  significantly  to  turbulence  activity  in  the  upper  troposphere  and  lower 
stratosphere,  at  altitudes  10-15  km,  and  several  accidents  in  high-altitude  aircraft  flights  have 
been  attributed  to  unexpected  mountain-wave-induced  stratospheric  turbulence  (Eckermann 
et  al.  [10]). 

From  a  theoretical  point  of  view,  the  generation  of  internal  gravity  waves  by  stratified  flow 
over  topography  is  a  classical  problem  of  fluid  dynamics.  Early  work  used  linear  models  to 
study  the  formation  of  two-dimensional  disturbances  on  the  lee  side  of  the  topography  -  the 
so-called  lee-wave  problem  -  under  perfectly  uniform  background  flow  conditions  (constant 
wind-velocity  and  buoyancy-frequency  profiles).  For  these  idealized  flow  conditions,  and  on 
the  further  assumption  of  steady  flow,  there  is  also  a  nonlinear  analytical  model  (Long  [17]) 
that  has  been  used  to  explore  wave  breaking  of  two-dimensional  disturbances  over  finite- 
amplitude  topography  (see  Baines  [1]  for  a  comprehensive  review). 

In  the  atmosphere,  however,  the  buoyancy  frequency  and  the  wind  speed  are  far  from  uni¬ 
form.  Typical  profiles,  from  balloon  measurements  in  the  lee  of  Mt.  Washington,  are  shown 
in  Fig.  1 .  The  mean  profiles  vary  significantly  with  altitude  near  the  tropopause,  the  border¬ 
line  between  the  troposphere  and  the  stratosphere.  In  this  instance,  the  tropopause  is  at  about 
12  km  from  the  ground  and  the  mean  buoyancy  frequency  increases  by  roughly  a  factor  of 
2  there.  In  addition,  the  buoyancy-frequency  profile  and  to  a  lesser  extent  the  wind  speed 
exhibit  appreciable  fluctuations  on  a  lengthscale  of  1-2  km  about  their  mean  values. 

In  more  recent  theoretical  work,  attempts  have  been  made  to  account  for  the  effect  of  non- 
uniform  background  buoyancy  frequency.  Extending  Long’s  model,  Durran  [9]  studied  non¬ 
linear  steady-state  hydrostatic  disturbances  over  finite-amplitude  topography  for  buoyancy- 
frequency  profiles  with  two  layers  of  uniform  stratification,  the  interface  modeling  the  tropopause. 
Depending  on  the  height  at  which  the  buoyancy-frequency  jump  is  positioned,  the  streamline 
deflections  can  be  maximized/minimized  over  the  topography,  thus  enhancing/suppressing 
wave  breaking  in  comparison  with  the  predictions  of  Long’s  model.  This  tuning/detuning 
effect  persists  for  continuous  buoyancy-frequency  profiles  with  a  smooth  transition  at  the 
tropopause,  as  demonstrated  by  Davis  [8]. 

For  the  purpose  of  forecasting  mountain  waves  in  the  field,  Eckermann  and  co-workers 
[7]  recently  devised  a  hybrid  methodology,  combining  linear  theory  to  estimate  the  near¬ 
field  response  close  to  the  topography,  with  a  ray-tracing  technique  that  accounts  for  gradual 
variations  of  the  buoyancy  frequency  and  wind  speed  in  the  far  field. 

In  the  present  work,  we  study  a  model  of  stratified  flow  over  topography  that  accounts 
for  the  tropopause  as  well  as  possible  oscillatory  variations  of  the  background  buoyancy  fre- 
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Fig.  1  Balloon  measurements  in  the  lee  of  Mt.  Washington  (a)  buoyancy  frequency,  and  (b)  flow  speed  as  a  function  of  altitude 
(courtesy  of  Dr.  R.  Beland). 


quency.  Thfe  potential  significance  of  the  latter  was  pointed  out  in  Prasad  and  Akylas  [20] 
using  an  asymptotic  theory  for  stratified  flow  with  constant  mean  buoyancy  frequency  near 
the  hydrostatic  limit.  Under  these  special  flow  conditions,  the  presence  of  small  sinusoidal 
variations  in  the  buoyancy  frequency  can  have  a  dramatic  effect  owing  to  a  resonance  mech¬ 
anism  brought  about  by  the  coupling  of  the  internal  gravity  waves  with  the  background  flow 
variations.  As  suggested  by  Phillips  [19],  this  wave-background-flow  interaction  causes  trap¬ 
ping  of  the  disturbance  close  to  the  topography;  as  a  result,  near  certain  critical  flow  speeds, 
the  response  behaves  similarly  to  resonant  flow  of  finite  depth,  exhibiting  significantly  in¬ 
creased  gravity-wave  activity  over  the  topography  as  well  as  upstream  influence  in  the  form 
of  solitary  waves  and  bores. 

Here  we  explore  this  resonance  mechanism  under  more  realistic  flow  conditions  by  solv¬ 
ing  the  Euler  equations  using  numerical  methods.  It  turns  out  that  periodic  variations  of  the 
buoyancy  frequency  continue  to  play  an  important  part  even  when  the  mean  buoyancy  profile 
is  not  constant  due  to  the  presence  of  the  tropopause,  although  strictly  the  theory  of  Prasad 
and  Akylas  [20]  is  no  longer  valid.  As  the  lengthscale  of  the  assumed  buoyancy-frequency 
oscillations  is  comparable  to  the  vertical  wavelength  of  the  induced  gravity  waves,  the  effects 
of  such  short-scale  variations  cannot  be  captured  by  ray  tracing. 

Preliminaries 

We  shall  use  the  same  non-dimensional  formulation  as  in  Prasad  and  Akylas  [20]  for  in- 
viscid  incompressible  flow  over  localized  topography  with  horizontal  lengthscale  L  and  peak 
amplitude  H.  Denoting  the  upstream  flow  speed  by  Uq  and  by  Nq  a  characteristic  value  of  the 
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Brunt-Vaisala  (buoyancy)  frequency,  the  flow  is  governed  by  three  independent  parameters 


Uq  HNq  NqU0 

^  =  777,  e  =  ~fT'  P  = 

NqL  Uq  g 


(1) 


g  being  the  gravitational  acceleration.  The  Boussinesq  parameter  /3  is  a  measure  of  strati¬ 
fication;  the  long-wave  parameter  p  and  the  amplitude  parameter  e  control  dispersive  and 
nonlinear  effects,  respectively. 

Employing  L  and  Uq/Nq  as  the  characteristic  lengthscales  along  the  horizontal  (x-)  and 
vertical  (y-)  directions,  respectively,  and  L/Uq  as  the  timescale,  the  governing  equations  (in 
the  Boussinesq  limit,  j3  — >0)  for  the  velocity  field  u  =  ( u ,  pv),  and  the  perturbations  of  density 
pp  and  pressure  p  from  hydrostatic  equilibrium,  are 


V-u  =  0,  (2) 

pt  +  u-  Vp—  N2v  =  0,  (3) 

11,  +  u-Vu  =  -(px,p~2(p  +  py)) .  (4) 

Here  N(y )  denotes  the  background  Brunt-Vaisala  frequency,  defined  by 

py  =  -ppN2,  (5) 


in  terms  of  the  background  density  p  (y) .  Moreover,  the  boundary  condition  on  the  topography 
y  =  eh(x)  reads 

v  =  £u-j-  ( y  =  eh ).  (6) 

ax 

In  the  present  study,  it  will  be  assumed  that  the  Brunt-Vaisala  frequency  consists  of  a 
mean  profile  with  oscillations  superposed  on  it.  Denoting  by  A  the  typical  lengthscale  of  the 
oscillations,  this  introduces  an  additional  parameter  which  we  take  to  be 


Upn 

AM)' 


(7) 


In  Prasad  and  Akylas  [20],  the  parameter  F  plays  the  role  of  a  Froude  number  as  it  defines 
the  critical-flow  regime  in  which  resonance  occurs,  in  analogy  with  resonant  flow  over  topog¬ 
raphy  in  a  channel  of  finite  depth  (Grimshaw  and  Smyth  [13],  Grimshaw  and  Yi  [12]). 

Specifically,  the  theory  of  Prasad  and  Akylas  [20]  applies  to  the  nearly  hydrostatic  re¬ 
sponse  (p  1)  in  the  case  that  small  sinusoidal  oscillations  of  wavelength  A  are  superposed 
on  a  uniform  mean  buoyancy  profile: 

M2(y)  =  1  +  jU2  (gisin(2Fy)  +  <??cos(2Fy)} .  (8) 


For  this  choice  of  N,  the  background  flow  is  nearly  uniformly  stratified  and,  according  to 
Long’s  model,  the  steady-state  response,  to  leading  order  in  p,  is  a  columnar  sinusoidal  dis¬ 
turbance  with  wavelength  equal  to  2k  in  the  vertical  direction;  the  weak  background-flow 
variations  in  (8),  small  deviations  from  the  hydrostatic  limit,  and  unsteady-evolution  effects 
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merely  modulate  this  columnar  disturbance.  Accordingly,  the  perturbation  streamfunction  \j/, 
in  terms  of  which  u  =  1  +  \jfy,  v  =  - y/x,  is  posed  as 

¥  =  {A(x,Y,Tyy  +  cc}  +  0(Li2),  (9) 

where  cc  denotes  the  complex  conjugate,  and  the  asymptotic  theory  furnishes  an  evolution 
equation  for  the  complex  envelope  A  of  the  columnar  disturbance  as  a  function  of  the  stream- 
wise  coordinate  x,  the  ‘stretched’  vertical  coordinate  Y  =  fi2y  and  the  ‘slow’  time  T  =  fi2t. 

From  Phillips  [19],  resonance  is  expected  to  occur  when  the  wavelength  of  the  background- 
flow  variations  is  about  one  half  of  the  response  wavelength,  so  F  &  1 .  Setting  then 

F  =  1  +  on2,  (10) 

a  =  0(  1)  being  a  resonance  detuning  parameter,  it  follows  that  the  linear  (e<C  1)  hydrostatic 
response  is  governed  by  the  evolution  equation 

AT-iAxY-\e2iaYQA*x  =  0,  (11) 

where  A*  is  the  complex  conjugate  of  A  and  Q  =  qi~  iq l ;  the  last  term  of  this  equation 
accounts  for  the  coupling  of  the  induced  columnar  disturbance  with  the  mean-flow  variations. 
In  addition,  from  (6)  and  (9),  the  linearized  boundary  condition  on  the  topography  is 

a  =  ~h{x)  (Y  —  0),  (12) 

where  a  is  the  real  part  of  A. 

The  critical  Froude  number  at  which  resonance  occurs  can  be  found  by  solving  equation 
(1 1)  subject  to  (27)  for  the  steady-state  linear  hydrostatic  response.  To  this  end,  we  write 

A=Ae,w  (13) 

so  that,  upon  substitution  into  (11),  A  satisfies  an  evolution  equation  with  constant  coeffi¬ 
cients.  The  steady  version  of  this  equation  corresponds  to  the  following  equation  system 

Oa  +  by -\qia  +  \q\b  =  0,  (14a) 

ab  -  aY  +  \q\a+  \qib  =  Q  (14b) 

in  terms  of  the  real  and  imaginary  parts  of  A  =  a  +  ib.  Moreover,  a  must  satisfy  the  boundary 
condition  (27).  The  solution  is  found  to  be 

a  =  -*Lh(x)e~sY,  b  =  ^h(x)  ^e~sY ,  (15) 

2  2  qi  +  4cr 

where 

5,2  =  —  cj2  +  T5  (?1  +?!)  •  (16) 

From  (15)  and  (16),  it  is  clear  that  the  steady-state  response  becomes  singular,  and  hence 
resonance  occurs,  when 

<r  =  -??2, 


q\  >  0. 


(17) 


(18) 


Combined  with  (10),  this  in  turn  determines  the  critical  Froude  number  F  —  Fcr jt, 

^crit  =  1  —  4?2/t  j  q\  >  0 


at  which  the  flow  is  resonant  according  to  the  linear  hydrostatic  theory. 

Prasad  and  Akylas  [20]  examined  the  implications  of  this  resonance  based  on  the  evolution 
equation  satisfied  by  the  complex  envelope  A  -  the  nonlinear,  weakly  dispersive  analogue  of 
(11)-  which  is  valid  asymptotically  when  small  sinusoidal  variations  are  present  on  an  other¬ 
wise  uniform  Brunt-Vaisala  frequency  in  nearly  hydrostatic  nonlinear  flow  over  topography] 
Here  we  explore  the  role  that  this  resonance  plays  under  less  restricted  flow  conditions  by 
solving  the  full  Euler  equations  (2)-(4)  numerically. 

Uniform  mean  buoyancy-frequency  profile 


We  have  constructed  a  numerical  model  for  solving  the  Euler  equations  (2)-(4)  for  strat¬ 
ified  flow  over  topography  with  uniform  wind  speed  but  general  buoyancy-frequency  profile 
far  upstream.  Since  the  flow  domain  is  unbounded,  a  major  issue  is  how  to  avoid  artificial  re¬ 
flections  at  the  boundaries  of  the  computational  domain.  For  this  purpose,  the  model  employs 
absorbing  viscous  layers  at  the  upper  and  lateral  boundaries.  A  description  of  the  numerical 
procedure  along  with  details  of  implementation  can  be  found  in  the  Appendix. 

In  all  our  computations,  we  have  used  the  algebraic  topography  profile,  also  known  as 
‘Witch  of  Agnesi’, 

(19) 


and  have  taken  £  =  0.6;  for  uniformly  stratified  (N=  1)  hydrostatic  flow  over  the  topography 
(19),  according  to  Long’s  model,  nonlinear  effects  come  into  play  roughly  above  this  value 
of  £  which,  however,  is  still  well  below  the  critical  value  of  £  =  0.85  required  for  overturning 
(wave  breaking). 

Our  first  choice  of  flow  conditions  mimics  those  considered  in  Prasad  and  Akylas  [20], 
namely  the  upstream  buoyancy-frequency  profile  is  given  by  (8)  with  fi2q\  =  0.25  cos0, 
Hzq2  =  0.25  sin0  (—it  <§  <k).  (The  sinusoidal  oscillations  are  terminated  at  y  =  12  it  and 
above  this  height  the  buoyancy  frequency  is  held  constant  at  N=  1;  the  numerical  stratifica¬ 
tion  profile  thus  comprises  12F  oscillation  periods,  in  contrast  to  the  asymptotic  theory  that 
assumes  infinitely  many  oscillations.)  According  to  (18),  for  this  choice  of  ji2q\  and  ji2q2, 
resonance  is  possible  when  —it/ 2  <  <j)  <  jt/2;  for  any  (j)  in  this  range,  the  critical  Froude 


number  then  is  given  by 


Fcrit 


sin  <j) 
16~’ 


(20) 


As  a  first  example,  we  take  0=0  and  F  =  1,  corresponding  to  exact  resonant  conditions 
according  to  (20).  Fig.  2(a)  shows  streamline  patterns  for  this  set  of  parameters  and  n  =  \/6 
at  t  =  50  and  t  =  150,  while  Fig.  2(b)  contrasts  the  response  at  t  —  150  against  the  solution 
we  would  have  obtained  at  the  same  time  in  the  absence  of  buoyancy-frequency  oscillations 
far  upstream;  by  t  =  150,  the  response  for  constant  N  has  essentially  reached  the  steady  state 
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predicted  by  Long’s  model.  It  is  evident  that  the  variations  of  the  Brunt-Vaisala  frequency 
alter  the  nature  of  the  flow  dramatically:  the  gravity-wave  activity  is  markedly  increased  over 
the  topography,  the  streamline  peak  amplitudes  associated  with  the  resonant  response  at  t  = 
150  being  nearly  twice  as  large  as  those  of  uniformly  stratified  flow.  In  the  resonant  response, 
moreover,  we  see  the  emergence  of  upstream-wave  propagation  along  with  a  disturbance 
of  opposite  sign  forming  on  the  downstream  side  of  the  topography.  This  behavior,  which 
is  characteristic  of  resonant  flow  in  a  channel  of  finite  depth  (Grimshaw  and  Smyth  [13], 
Grimshaw  and  Yi  [12]),  was  also  found  in  Prasad  and  Akylas  [20]  and  must  be  attributed  to 
the  trapping  effect  caused  by  the  interaction  of  the  induced  disturbance  with  the  background 
buoyancy-frequency  oscillations.  In  the  asymptotic  theory,  the  response  was  tracked  in  terms 
of  the  scaled  time  T  =  fl2t,  corresponding  to  significantly  larger  values  of  t  than  those  in  Fig. 
2;  at  such  later  times,  typically,  wave  breaking  occurs  and  the  upstream  disturbance  evolves 
into  solitary  waves  or  bores. 

Fig.  3  shows  two  snapshots  of  the  response  for  fi  —  0.05,  F  =  1  and  (j>  =  n/3  instead  of 
<j)  =  0.  From  (18),  this  change  of  the  phase  <p  lowers  the  value  of  the  critical  Froude  number 
to  Fcrjt  =  0.957  so  F  =  1  is  supercritical  under  the  present  flow  conditions.  As  a  result,  the  re¬ 
sponse  in  Fig.  3  is  quite  different  from  that  shown  in  Fig.  2(a):  the  streamline  peak  amplitudes 
are  generally  smaller,  there  is  no  upstream-wave  propagation,  and  apparently  steady  state  is 
approached  as  the  transients  are  swept  downstream.  Fig.  3  also  shows,  for  comparison,  the 
hydrostatic  steady-state  response  corresponding  to  the  present  flow  conditions,  as  obtained 
following  the  numerical  procedure  of  Davis  [8];  this  confirms  that  the  unsteady  response  in¬ 
deed  reaches  steady  state.  On  the  other  hand,  by  lowering  the  Froude  number  to  the  critical 
value  of  f'crit  —  0  .957  given  by  (18),  the  response  again  becomes  resonant  (see  Fig.  4).  (We 
remark  that  the  smaller  value  of  fi  =  0.05,  which  is  closer  to  the  hydrostatic  limit,  rather  than 
/i  =  1  /6,  was  used  in  Fig.  3,  only  to  facilitate  comparison  against  the  hydrostatic  steady  state; 
the  difference  between  the  resonant  and  non-resonant  behavior  is  controlled  by  F  and  is  not 
affected  by  this  choice  of  ju.) 

The  results  of  fully  numerical  simulations  reported  above  support  the  conclusions  reached 
by  the  asymptotic  theory  of  Prasad  and  Akylas  [20].  Superposing  sinusoidal  variations  on 
a  uniform  mean  buoyancy-frequency  profile  indeed  causes  trapping  of  gravity  wave  distur¬ 
bances  in  the  vertical  direction,  and  the  flow  behaves  as  if  it  were  in  a  channel  of  finite 
depth.  The  response  turns  out  to  be  very  sensitive  to  small  changes  in  the  background-flow 
conditions  and,  near  the  critical  Froude  number,  features  increased  gravity-wave  activity  and 
upstream  influence.  In  the  following  section,  we  explore  the  role  of  this  resonance  mechanism 
when  the  buoyancy  profile  takes  into  account  the  tropopause. 

Effects  of  the  tropopause 

We  shall  model  the  transition  from  the  troposphere  to  the  stratosphere  by  the  dimension¬ 
less  mean  buoyancy-frequency  profile 


(N)  —  1.5  +  0.5tanh2(y  —  d), 


(21) 
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(b) 

Fig.  2  Resonant  response  for  F  =  1 ,  <p  =  0  and  p  =  1/6  (a)  evolution  in  time:  dashed  line  t  =  50,  solid  line  t  =  1 50;  (b)  comparison 
between  the  resonant  response  (solid  line)  and  the  solution  for  uniform  buoyancy  frequency  profile  (dashed  line)  at  t=  150. 


where  d  is  a  parameter  that  controls  the  tropopause  height.  Consistent  with  Fig.  1(a),  the 
profile  (56)  is  such  that  (N)  varies  smoothly  from  1  in  the  troposphere  (y  d)  to  2  in  the 
stratosphere  (y^>d).  Referring  to  Fig.  1,  a  typical  value  of  the  mean  buoyancy  frequency 
in  the  troposphere  is  Nq  =  0.01  sec-1  and  the  wind  speed  is  in  the  range  5-30  m/sec,  so 
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Fig.  3  Non-resonant  response  for  F=l,  0  =  7t/3  and  /i=0.05  at/ = 50  (dotted  line)  and /=  100  (dashed  line).  The  corresponding 
hydrostatic  steady-state  solution  (solid  line)  is  also  shown  for  comparison. 


Fig.  4  Time  evolution  of  resonant  response  for  F =0.957,  <j>  =  n/ 3  and  n  =  1/6:  dashed  line  t  =  50,  solid  line  t  —  1 50. 


the  characteristic  lengthscale  Uq/Nq  ~0.5  —  3  km.  Generally,  the  atmospheric  tropopause 
height  varies  from  10  to  14  km  which,  in  terms  of  the  dimensionless  parameter  d,  would  then 
correspond  to  d  r sj  5-25. 
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The  steady-state  hydrostatic  response  for  a  two-layer  buoyancy-frequency  profile  with 
a  finite  jump  at  the  interface  of  the  two  layers  was  studied  by  Durran  [9]  and  later  Davis 
[8]  considered  the  same  problem  for  a  continuously  varying  buoyancy  frequency  as  in  (56). 
They  both  find  specific  tropopause  heights  for  which  the  hydrostatic  steady-state  response 
exhibits  enhanced  or  reduced  gravity-wave  activity.  For  e  =  0.6,  for  example,  such  ‘tuned’ 
and  ‘detuned’  gravity-wave  flow  fields  occur  when  the  tropopause  is  placed  in  the  vicinity 
of  J  =  21  and  d  —  25,  respectively.  In  fact,  consistent  with  Durran  [9]  and  Davis  [8],  we  find 
multiple  ‘tuned’  and  ‘detuned’  responses  that  recur  periodically  with  changing  d,  the  lowest 
value  of  d  that  corresponds  to  a  ‘tuned’  response  being  in  the  neighborhood  of  d — 6. 

For  two-layer  buoyancy-frequency  profiles  with  sharp  interface,  Durran  [9]  and  Davis  [8] 
also  noted  that  their  numerical  procedures  failed  to  converge  to  a  hydrostatic  steady-state 
solution  at  certain  tropopause  heights.  While  Durran  [9]  suggested  that  this  anomaly  could 
indicate  a  nonlinear-resonance  phenomenon,  Davis  [8]  attributed  the  lack  of  convergence  to 
numerical  difficulties. 

For  the  purpose  of  validating  and  clarifying  the  steady-state  results  of  Durran  [9]  and 
Davis  [8],  we  have  carried  out  unsteady  numerical  simulations  using  the  buoyancy-frequency 
profile  (56)  close  to  the  hydrostatic  limit  (/i=0.05).  We  find  that  the  unsteady  response  gener¬ 
ally  approaches  the  steady  state  predicted  by  the  previous  studies;  a  typical  example  is  shown 
in  Fig.  5  for  d  —  3.  We  also  investigated  the  nature  of  the  response  for  buoyancy  profiles  with 
sharp  transition  at  the  tropopause  for  the  values  of  d  at  which  Durran  [9]  and  Davis  [8]  ex¬ 
perienced  difficulties  converging  to  a  hydrostatic  steady  state.  Near  these  special  values  of  d 
(d=4.5, 11, 17.3,23.5,...),  it  turns  out  that  the  steady-state  response  approaches  the  uniform 
stream  far  upstream  and  downstream  much  more  slowly  than  Long’s  solution,  and  a  suffi¬ 
ciently  large  computational  domain  is  required  along  x  in  order  to  achieve  convergence  (see 
Fig.  6).  Also  our  numerical  simulations  indicate  that  the  corresponding  unsteady  response 
tends  to  steady  state,  albeit  relatively  slowly,  and  there  is  no  sign  of  a  nonlinear  resonance. 

Next,  we  superpose  on  the  mean  buoyancy  profile  (56)  sinusoidal  oscillations  similar  to 
those  used  in  (8)  for  the  uniform  mean  profile: 

N2(y)  =  {1.5  +  0.5  tanh  [2(y  —  d)]}2  +  /i2  {q\  sin(2Fy)  +  g2c°s(2Fy)}  (22) 

with  ji2q\  —  0.25  cos <j),  H2qi  =  0.25  sin0  as  before.  Fig.  7  shows  the  unsteady  response  for 
H  —  1/6,  0  =  0  and  F  —  1  at  t  —  50  and  t  =  150  when  the  tropopause  is  placed  at  d  —  21. 
Note  that,  for  uniform  mean  buoyancy  frequency  equal  to  1,  this  choice  of  flow  parame¬ 
ters  would  correspond  to  resonant  conditions  in  the  troposphere  (y<^d)  according  to  (20), 
but  for  the  buoyancy-frequency  profile  (22),  which  includes  the  tropopause,  strictly  the  reso¬ 
nance  equation  (20)  is  no  longer  valid.  Nevertheless,  upon  comparing  Fig.  7  with  Fig.  2,  the 
salient  features  of  resonant  response  identified  earlier  -  increased  wave  amplitudes  over  the 
topography  and  upstream-wave  formation  -  are  still  present,  suggesting  that  the  resonance 
mechanism  persists  when  the  tropopause  is  taken  into  account;  even  when  trapping  occurs 
in  the  troposphere  only,  resonance  is  still  possible.  Finally,  Fig.  8  shows  the  response  for  the 
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Fig.  5  Comparison  of  the  response  for  tropopause  height  d  =  3  and  ji  =  0.05  at  t  —  50  (dashed  line)  with  the  corresponding 
hydrostatic  steady-state  solution  (solid  line).  At  this  time,  steady  state  has  essentially  been  reached  below  y= 3. 


x 


Fig.  6  Contrast  between  Long’s  solution  (solid  line)  and  the  steady-state  hydrostatic  response  for  d  =  23.25  computed  using  a 
domain  of  width  x=  ±50  (dotted  line)  and  *=  ±100  (dashed  line). 


same  flow  parameters  as  in  Fig.  7  but  with  <p  —  n/3.  As  expected,  this  response  is  not  resonant 
and  approaches  steady  state  as  the  non-resonant  case  shown  in  Fig.  3. 
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Fig.  7  Time  evolution  of  resonant  response  for  F  =  1 ,  ^  =  0  and  /t  =  1/6  in  the  presence  of  a  tropopause  at  d  =  2 1:  dashed  line 
t = 50,  solid  line  t  =  150. 


x 

Fig.  8  Response  at  t  =  150  for  the  same  flow  parameters  as  in  Fig.  7  but  with  b  =  n/3. 


The  results  from  numerical  simulations  of  stratified  flow  over  topography  reported  here 
suggest  that  periodic  oscillations  in  the  background  buoyancy  frequency  cannot  be  ignored. 
A  study  of  more  general  (quasi-periodic,  unsteady)  short-scale  variations  in  the  background 
flow  conditions  is  currently  under  way. 
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Numerical  method  for  unsteady  computations 

The  unsteady  Euler  equations  (2)-(4)  were  solved  using  the  second-order  projection  method 
developed  by  Bell  and  co-workers  [3-5].  Lamb  [15]  implemented  this  procedure  to  uniformly 
stratified  finite-depth  flow  over  locally  confined  topography  and  Skopovi  [21]  allowed  for 
variations  in  the  buoyancy-frequency  profile.  Here,  we  apply  the  numerical  scheme  to  strat¬ 
ified  flow  of  large  depth  over  topography  by  introducing  viscous  layers  that  dissipate  dis¬ 
turbances  leaving  the  inviscid  region  and  thus  avoid  reflections  from  the  boundaries  of  the 
computational  domain. 

For  this  purpose,  we  solve  the  governing  equations  (2)— (4)  with  the  viscous  term  v(ii2yixx + 
Ujy)  added  to  the  right-hand  side  of  the  momentum  equation  (4).  Here  v  denotes  a  spatially 
varying  dimensionless  viscosity  that  is  distributed  within  viscous  layers  around  the  inviscid 
region.  Specifically,  near  the  upper  boundary, 

v  =  (23) 

so  v  increases  gradually  from  zero  at  the  lower  edge  (y  =  yo)  to  v  =  0. 1  ji~]  at  the  upper 
edge  (y  =  60)  of  the  viscous  layer.  The  upper  edge  coincides  with  the  top  boundary  of  the 
computational  domain.  In  the  side  layers,  to  ensure  numerical  stability  at  the  inflow  and 
outflow  boundaries,  the  viscosity  is  first  increased  in  a  way  analogous  to  (55),  then  is  held 
constant  at  O.ljU-1  for  M  number  of  grid  points  and  finally  is  brought  back  to  zero  again 
sinusoidally  over  the  span  of  M  points  (see  Skopovi  [22]  for  details). 

In  all  simulations,  we  used  yo  =  41,  while  M= 50  for  computations  with  fi  =  1/6  and  M= 
25  for  runs  with  ju  =  0.05.  A  typical  contour  plot  of  the  viscosity  profile  is  presented  in  Fig. 

1 1  from  which  it  is  evident  that  nearly  two  thirds  of  the  computational  domain  are  occupied 
by  absorbing  layers.  In  order  to  accommodate  the  bottom  boundary  of  the  computational 
domain  where  we  enforce  the  inviscid  boundary  condition  (i.e.  zero  velocity  perpendicular  to 
the  topography),  viscosity  in  the  lateral  absorbing  layers  is  linearly  increased  in  the  vertical 
direction  from  zero  at  the  topography  to  the  layer  value  over  20  grid  points,  as  indicated  by 
curving  contours  in  the  bottom  left  and  right  comers  of  Fig.  11. 

The  numerical  method  of  solution  was  tested  in  various  ways.  First  we  verified  that  the 
nearly  hydrostatic  unsteady  response  for  uniformly  stratified  flow  over  topography  converges 
towards  Long’s  steady-state  solution.  This  was  a  rigorous  test  of  the  performance  of  the  upper 
viscous  layer,  as  in  the  absence  of  oscillations  in  the  stratification  profile,  disturbances  are  not 
trapped  so  all  energy  is  radiated  upwards.  Secondly,  as  a  test  of  the  lateral  absorbing  layers, 
we  compared  the  response  using  viscous  layers  at  x  =  ±18  against  that  obtained  when  the 
layers  are  positioned  at  *=±10.  Typical  results  are  shown  in  Fig.  10  from  which  it  is  apparent 
that  the  layers  are  effectively  absorbing  waves  exiting  the  inviscid  region  without  impacting 
the  flow  field  within  the  region  of  interest.  We  remark  that  this  comparison  is  made  at  t=  150, 
which  corresponds  to  the  largest  time  at  which  we  present  our  unsteady  results.  Lastly,  we 
verified  that  our  numerical  results  did  not  vary  appreciably  when  the  spatial  resolution  was 
changed. 
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Fig.  10  Comparison  of  the  resonant  response  for  F=  1,  fi  =  1/6,  tp  =  0  and  t  =  150  when  absorbing  layers  are  positioned  at 
x=±10  (dashed  line)  with  the  response  obtained  for  the  same  set  of  parameters  when  the  inviscid  region  extends  to  x=  ±18 
(solid  line). 


Following  Lamb  [15],  we  have  used  a  quadrilateral  grid  with  the  higher  resolution  in  the 
vicinity  of  the  obstacle.  The  grid  is  linearly  stretched  along  both  the  horizontal  and  the  vertical 
such  that  the  height  of  the  bottom  cells  is  15  times  smaller  than  the  height  of  the  cells  at  the 
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upper  boundary,  while  the  width  of  the  middle  cells  is  15  times  smaller  than  the  width  of 
the  left-  and  right-boundary  cells  for  simulations  where  [i  =  1/6,  and  10  times  smaller  for 
runs  where  ji =0.05.  All  simulations  are  carried  out  on  a  numerical  domain  that  contains  245 
cells  in  the  vertical  direction.  In  the  horizontal  direction,  1200  cells  were  used  for  flows  with 
[1  =  1/6  and  800  cells  for  flows  with  [1=0.05.  In  all  runs,  the  width  of  the  numerical  domain 
is  x  £  [-30,30]  when  [i  =  1/6  and  x  £  [—15, 15]  when  [1  =  0.05.  The  only  exception  is  the 
simulation  in  Fig.  10  with  the  viscous  layers  starting  at  x=  ±10,  which  is  produced  with  860 
cells  in  the  horizontal  direction  and  the  lateral  boundaries  of  the  computational  domain  at 
x=±20.  Lastly,  in  regard  to  temporal  discretization,  upon  impulsively  accelerating  the  flow 
from  rest,  we  use  the  time  step  At =0.015  for  runs  with  [i  =  1/6  and  At=  0.02 with  [ i =0.05. 
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B.  Effects  of  Temporal  Variations  in  Wind  Velocity  on  the  Generation  of  Mountain 
Gravity  Waves 

The  vast  majority  of  previously  conducted  studies  on  the  subject  of  stratified  flows  over  to¬ 
pography,  primarily  for  mathematical  tractability  and  computational  cost,  assumes  a  constant 
background  velocity.  It  is  generally  known,  however,  that  a  typical  atmospheric  free  stream, 
in  addition  to  exhibiting  spatial  variations,  regularly  fluctuates  with  time.  In  this  study,  we  are 
interested  in  the  effect  of  these  temporal  modulations. 

On  the  theoretical  side,  one  of  the  first  comprehensive  efforts  to  account  for  time  variations 
in  the  wind  velocity  is  the  work  of  Bell  [6]  who  in  1975  derived  the  linear  steady-in-the-mean 
solution  for  a  horizontally  incident  and  sinusoidally  oscillating  background  velocity  with 
zero  mean.  He  established  that  the  equation  governing  his  model  is  nonlinear  with  respect 
to  the  basic  flow.  As  a  result,  internal  gravity  waves  in  his  model  were  produced  not  only  at 
the  fundamental  frequency,  but  also  at  all  of  its  harmonics.  Later,  Bannon  and  Zehnder  [2] 
constructed  a  model  that  is  more  appropriate  for  atmospheric  conditions  by  adding  a  steady 
component  to  the  sinusoidally  varying  free  stream  of  Bell  [6].  In  spite  of  the  fact  that  it  is 
restricted  to  a  hydrostatic  mountain,  their  linear  steady-in-the-mean  solution  encompassed  the 
Coriolis  acceleration  and  provided  perhaps  the  first  insight  into  the  importance  of  temporal 
modulations  in  the  atmospheric  context  by  indicating  that  the  instantaneous  mountain  drag 
may  be  larger  than  that  exerted  by  the  steady  wind.  We  also  cite  the  work  of  Hines  [14]  who 
allowed  the  steady  portion  of  the  background  velocity  to  vary  with  height,  but  slowly  so  that 
the  WKB  approximation  can  be  used  to  analyze  the  resulting  flow  field.  Although  Hines’s 
linear  analysis  is  the  first  substantial  effort  to  consider  simultaneously  shear  and  unsteadiness 
in  the  basic  flow,  it  is  rather  limited  to  a  monochromatic  mountain  profile. 

Among  computational  studies,  we  point  out  the  work  of  Lott  and  Teitelbaum  [18]  who 
investigated  the  formation  of  unsteady  linear  gravity  waves  in  a  free  stream  that  starts  from 
zero  and  returns  to  zero  after  a  finite  time.  Their  analysis  is  limited  to  a  monochromatic 
topography,  but  it  includes  effects  of  shear.  We  also  acknowledge  the  work  of  Eckermann 
and  co-workers  [7]  who,  for  the  purpose  of  forecasting  mountain  waves  in  the  field,  recently 
devised  a  hybrid  methodology.  This  approach  combines  linear  theory  to  estimate  the  near¬ 
field  response  close  to  the  topography  with  a  ray-tracing  technique  that  accounts  for  gradual 
variations  of  the  buoyancy  frequency  and  wind  speed  in  the  far  field. 

These  previously  conducted  studies  that  account  for  temporal  changes  in  the  background 
velocity  possess  several  limitations.  Primarily,  they  are  linear  in  nature  and  they  neglect  short- 
scale  buoyancy  frequency  variations.  In  the  work  of  Prasad  and  Akylas  [20]  and  later  Skopovi 
and  Akylas  [23],  these  mechanisms  have  been  found  to  play  a  substantial  role  in  shaping  the 
character  of  the  response  with  the  steady  free  stream.  At  this  point,  therefore,  a  natural  ques¬ 
tion  to  ask  is  how  does  the  flow  behave  in  the  absence  of  these  constraints.  In  an  attempt 
to  answer  this  question,  we  first  derive  the  linear  hydrostatic  steady-in-the-mean  solution  for 
sinusoidally  oscillating  background  velocity,  which  is  valid  for  an  arbitrary  Brunt-Vaisala 
frequency.  The  solution  and  details  of  its  numerical  implementation  are  discussed  in  the  fol¬ 
lowing  section.  Furthermore,  we  construct  a  nonlinear  computational  model  capable  of  sim- 
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ulating  the  flow  field  generated  by  the  time-varying  background  flow.  Implementation  details 
of  the  model  are  presented.  Equipped  with  these  two  tools,  we  explore  the  significance  of 
nonlinearity  by  examining  the  response  due  to  a  gradually  increasing  background  flow  from 
zero  to  a  constant  value  in  the  presence  of  uniform  stratification.  Following  this  effort,  we 
turn  our  attention  to  1-2  km  modulations  in  the  buoyancy  frequency.  We  explore  their  im¬ 
portance  by  considering  the  transient  response  to  the  basic  velocity  that  is  the  superposition 
of  a  uniform  stream  and  harmonic  oscillations.  Specifically,  we  concentrate  on  frequencies 
that  are  of  the  order  of  the  atmospheric  Brunt-Vaisala  frequency  and  finally  we  analyze  the 
behavior  at  periods  that  are  of  the  order  of  hours. 

Linear  Theory  for  Monochromatic  Variations 


We  shall  use  the  same  non-dimensional  formulation  as  in  Prasad  and  Akylas  [20]  for 
inviscid  incompressible  flow  over  localized  topography  with  horizontal  lengthscale  L  and 
peak  amplitude  H.  Denoting  the  upstream  flow  speed  by  C/q  and  by  No  the  characteristic 
value  of  the  Brunt-Vaisala  (buoyancy)  frequency,  the  flow  is  governed  by  three  independent 
parameters 


Uo  „_HN0  o_N0U0 

P  “  xt  r  ?  ^  jr r  J  P  ? 
NqL  Uo  g 


(24) 


g  being  the  gravitational  acceleration.  The  Boussinesq  parameter  ft  is  a  measure  of  strati¬ 
fication;  the  long-wave  parameter  j u  and  the  amplitude  parameter  e  control  dispersive  and 
nonlinear  effects,  respectively. 

Employing  L  and  Uq /Nq  as  the  characteristic  lengthscales  along  the  horizontal  (x-)  and 
vertical  ( z -)  directions,  respectively,  and  L/Uq  as  the  timescale,  the  governing  equations  (in 
the  Boussinesq  limit,  j3  ->0)  become 


( Tt+UlTx )  Ozz  +  MV**)  -rx  =  eJ(y,yzz  +  ^i2yxx),  (25a) 

(§'t+U'L)r+N2Vx  =  eJ{W^  (25b) 

where  J(a,b)  =  axbz  —  azbx.  In  these  expressions,  y  is  the  perturbation  stream  function, 
which  physically  denotes  the  deflection  of  a  streamline  from  the  hydrostatic  value  it  possesses 
far  upstream  of  the  obstacle.  It  is  physically  defined  as  ( u ,  fiw)  =  (U +eyz,  —£yx)  where  u  and 
w  are  respectively  horizontal  and  vertical  velocity  components.  Moreover,  U(t)  denotes  the 
background  flow  velocity,  N  represents  the  buoyancy  frequency  while  ej3r  is  the  pretrubation 
denisty. 

In  order  to  complete  the  mathematical  formulation  of  the  problem  under  consideration,  we 
need  to  specify  boundary  conditions.  Due  to  the  fact  that  the  topography  is  presumed  to  be 
locally  confined  in  the  horizontal  direction  and  the  fluid  is  considered  to  be  infinitely  deep, 
we  require 


y  =  0  and  r  =  0  at  x  =  ±°°  and  z  —  +°°. 


(26) 
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Along  the  bottom  boundary,  on  the  other  hand,  the  velocity  component  perpendicular  to  the 
barrier  is  zero,  which  in  turn  implies  that  the  scalar  product  of  the  velocity  and  the  unit  vector 
orthogonal  to  the  mountain  is  zero.  Algebraically,  this  can  be  written  in  terms  of  yf  as 

dh 

-V6c  =  ^-(f/  +  e%)  at  z  =  £h(x),  (27) 

where  h(x)  designates  the  topography  profile.  At  this  point,  it  is  simportant  to  emphasize  that 
in  all  simulations  in  this  study,  we  use  the  topography  profile 


T+^2 

popularly  known  as  the  algebraic  mountain  or  Witch  of  Agnesi. 

We  commence  the  derivation  of  the  linear  steady -in-the-mean  solution  that  accommodates 
any  stratification  by  defining  the  Fourier  transform  and  its  inverse  with  respect  to  x  as 


=£ 


y/e  lkxdx  and  $ 


-  r°°yreikjcdk. 


Taking  the  transform  (29)  of  our  governing  equations  (25)  in  the  limit  e  ->  0  then  yields 


^  +  (-k2!l2yr+y/zz)-N2k2yr  =  0. 


Following  Bell  [6],  (30)  can  be  mathematically  simplified  by  introducing  the  horizontal  co¬ 
ordinate  (E,)  traveling  with  the  time- varying  free  stream.  In  terms  of  x  and  U,  E,  is  defined 


E,  =  x-  f  U(x)dx. 
Jo 


The  relationship  between  yr  and  the  Fourier  transform  with  respect  to  E, ,  denoted  as  yf,  is  then 

yf=yfeik^U^dx.  (32) 

Substituting  (32)  into  (30)  gives  the  linearized  governing  equation  (25)  in  terms  of  y 'r  as 

Pa tt  ~  Vtt  -  N2k 2  yf  =  0.  (33) 

We  construct  the  bottom  boundary  condition  by  inserting  Fourier  transform  of  (27)  into 
(32)  to  obtain 

V=~JkJt  (/*/<! U(T)dZ)  at  z  =  0.  (34) 

Here,  h  designates  the  Fourier  transform  of  the  ridge.  Following  Bannon  and  Zehnder  [2],  at 
this  point,  we  restrict  our  considerations  to  the  background  velocity  of  the  form 

U  =  1  -t-Acos(o)oO  (35) 
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where  both  A  and  coo  are  constants.  With  this  basic  flow,  the  time  derivative  in  (34)  may  be 
eliminated  by  introducing  the  Bessel  functions  of  the  first  kind  (/„)  defined  by 


exp 


y 

2 


=  I  cfJniy) 


tl—~  ®o 


where  n  is  an  integer.  For  ln(g)  =  itp,  it  follows  that 


(36) 


giy sin <p  _  £  e^nJn{y). 

n=—oo 


(37) 


Substitution  of  (35)  and  (37)  into  (34)  with  y  =  kA /(Do  and  0  =  (Dot  leads  to  the  bottom 
boundary  condition 

f=-h  V  n(—\e at  Z  =  0.  (38) 

ntL~  k  \ODoJ 


Based  on  (38),  we  seek  the  solution  of  (33)  in  the  form 


+° 


Mk+(Oon)t 


(39) 


where  W  is  the  function  dependent  on  z,  k  and  n.  Substitution  of  (39)  into  (33),  reveals  that 
W  must  satisfy  the  second  order  ordinary  differential  equation  with  respect  to  z 


W'+fc2 


N2 


(k  +  (Qon)2 


•r 


w  =  o. 


(40) 


Before  discussing  a  varying  stratification,  we  consider  the  case  of  constant  Brunt-Vaisala 
frequency.  In  this  scenario,  the  solution  of  (40)  can  be  determined  analytically  and  it  takes 
the  form 


W  =  C\eimz  +  C2  e~imz 


(41) 


with  m  given  by 

_  [  k2N2 
m  (k  +  (Don)2 

In  (41),  Cj  and  C2  are  functions  of  k  and  n.  They  are  determined  from  (38)  and  the  top 
boundary  condition. 

The  form  of  the  upper  boundary  condition  for  each  k  and  n  depends  on  whether  (41)  is 
exponential  or  oscillatory  in  nature.  It  is  exponential  if  N2/ (k+  (Don)2  -  jl2  <  0.  In  this  case, 
the  portion  of  (41)  that  grows  with  height  is  physically  unrealistic  and  it  is  neglected  by 
setting  €2  =  0.  The  solution  of  (41)  then  becomes 


kV 


(42) 


W  =  Ceimz  with 


m  ■  i 


k 

k  +  (Don 


yjli2(k  +  (0on)2  -N2 


(43) 
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with  C=Ci .  If  N2/ (k  +  (Don)2  —  fi2  >  0,  (41)  is  oscillatory.  In  this  regime,  (39)  indicates  that 
is  the  superposition  of  plane  waves  with  the  frequency  co=—(k+  (Don)  and  the  vertical 
wavenumber  that  is  either  a  positive  or  negative  quantity  (42).  The  proper  sign  of  m  for  each 
k  and  n  is  determined  from  the  radiation  condition,  which  states  that  d(D/dm  >  0.  In  terms  of 
(Dq,  k  and  n,  this  can  be  expressed  as 


d(D  _  m{k  +  (Oon) 
dm  m2  +  [l2k 2  > 


(44) 


Due  to  the  fact  that  the  denominator  of  this  inequality  is  always  greater  than  zero,  the  group 
velocity  is  positive  and  and  all  energy  is  radiated  outwards  if  the  sign  of  m  coincides  with  that 
of  k+  (Don.  Consequently,  in  the  oscillatory  regime,  (41)  may  be  written  as 


W  =  Ceimz  with 


m  — J — — J N2  —  ll2(k-\- (Don)2 

k  +  (DonV  K  ' 


(45) 


where  C  is  either  Ci  or  C2  depending  on  the  sign  of  k+  (Don.  In  (43)  and  (45),  C  is  determined 
from  (38)  to  be 


C  =  - 


k  +  (Don 
k 


(46) 


The  linear  steady-in-the-mean  solution  for  constant  N  and  U  given  by  (35)  may  now  be 
written  as 


ei[kx+mz-kA/ Wq  si^Oo^+OJont]^ 


(47) 


where  W  is  determined  from  (43)  and  (45). 

It  can  be  easily  verified  that  for  the  uniform  basic  velocity  ( A  =0)  and  N=  1,  the  hy¬ 
drostatic  limit  (ju  ->0)  of  (47)  reduces  to  the  linear  solution  of  Long  [17].  Moreover,  in  the 
absence  of  the  time-invariant  component  of  the  free  stream  (35),  (47)  is  equivalent  to  the 
equation  governing  the  model  of  Bell  [6].  Lastly,  (47)  is  analogous  to  the  solution  of  Bannon 
and  Zehnder  [2]  when  the  Coriolis  acceleration  in  their  model  is  neglected. 

The  solution  (47)  can  be  computed  by  using  Gaussian  quadrature  to  evaluate  the  integral 
for  each  n.  In  that  regard,  it  is  important  to  note  that  k =0  is  an  integrable  singularity.  In  fact, 
using  FHopital’s  rule,  one  can  demonstrate  that  the  limit  0  of  the  integrand  is  finite  for 
all  n  regardless  of  stratification. 

The  solution  of  (40)  for  a  varying  buoyancy  frequency  must  be  determined  numerically. 
While  the  bottom  boundary  condition  is  given  by  (38)  regardless  of  the  stratification,  the 
upper  boundary  condition  for  constant  N  can  still  be  used  if,  following  Davis  [8],  we  assume 
that  above  some  height  z=z^,,  N=Nca.  For  numerical  simplicity,  we  additionally  restrict  our 
considerations  to  the  hydrostatic  response  (ju  — >0).  With  these  two  assumptions,  the  upper 
boundary  condition  for  constant  N  is  given  by 


w  =  , 


(48) 
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and  its  ^-derivative  is 


WZ  =  C 


i\k\N„ 

k+a>on 


o‘t+a Qnz°° 


(49) 


We  solve  (40)  by  first  guessing  the  value  of  C  in  (48)  and  (49).  From  the  knowledge  of  W 
and  Wz  at  z=Zoo,  we  then  use  the  fourth  order  Runge-Kutta  method  to  attain  W  at  z=0.  The 
difference  between  this  computed  value,  denoted  here  as  Wo,  and  that  required  by  the  bottom 
boundary  condition  (38)  gives  the  error  associated  with  our  guess,  namely, 


£err(C)=W0(C)  +  ^p. 


(50) 


We  update  this  guess  via  Newton-Raphson  scheme.  In  particular,  the  new  guess  is  given  by 


Cnew  =  C- 


Eerr  (Q 
E'errW 


(51) 


Here,  E'err  is  the  approximation  of  the  error  function’s  derivative  at  C  obtained  from 


E'eAC)  = 


Eerr(C  +  S)-Eerr{C) 


(52) 


where  5  is  a  small  perturbation  to  C.  Once  we  obtain  W  for  each  k  and  n,  we  evaluate  (47) 
using  Gaussian  quadrature  analogous  to  the  case  of  constant  N. 

In  all  simulations  presented  in  this  document,  we  use  A:  €  [—15, 15]  with  the  grid  spacing 
of  0.025.  Moreover,  we  choose  the  range  of  n  in  such  a  way  that  for  all  values  of  n  outside 
this  range  Jn(kA/(Oo )  <  10-5.  In  all  computations  involving  a  varying  stratification,  we  use 
8  =  10“5  and  the  initial  guess  C  =  —  (k  +  0)Qn)/k.  We  iterate  the  Newton-Raphson  scheme 
until  \Eerr\  <  10-3. 

The  numerical  method  for  non-uniform  N  was  tested  in  two  ways.  First,  we  verified  that 
for  N  =  1  and  A  =0,  the  solution  coincides  with  that  of  Long  [17].  Moreover,  for  A—0  and 
a  two-layer  buoyancy  frequency,  we  confirmed  the  agreement  with  the  linear  (£-»0)  limit  of 
the  response  of  Durran  [9]. 

Numerical  method  for  unsteady  computations 

In  this  study,  we  solve  the  unsteady  Euler  equations  (25)  using  the  second-order  projection 
method  developed  by  Bell  and  co-workers  [3-5].  Lamb  [16]  implemented  this  procedure  to 
uniformly  stratified  finite-depth  flow  over  locally  confined  topography  and  Skopovi  [21]  al¬ 
lowed  for  variations  in  the  buoyancy-frequency  profile.  In  the  further  development,  Skopovi 
and  Akylas  [23]  applied  the  numerical  scheme  to  stratified  flow  of  large  depth  over  topogra¬ 
phy  by  introducing  viscous  layers  that  dissipate  disturbances  leaving  the  inviscid  region  and 
thus  avoid  reflections  from  the  boundaries  of  the  computational  domain.  Here,  we  addition¬ 
ally  extend  the  scheme  to  encompass  temporal  variations  in  the  background  velocity. 

We  implement  the  time  modulations  of  the  free  stream  velocity  U  by  adding  the  forcing 
function  b(t )  to  the  right-hand  side  of  the  horizontal  component  of  the  momentum  equation. 
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This  term  physically  denotes  the  local  time  rate  of  change  of  the  spatially  uniform  horizontal 
background  velocity.  Therefore,  accommodating  a  particular  temporal  variation  in  the  basic 
flow,  from  the  numerical  standpoint,  merely  entails  specifyimg  this  forcing  function.  In  this 
study,  we  are  concerned  with  a  free  stream,  U,  that  is  the  sum  of  a  homogeneous  profile  and 
a  monochromatic  time- varying  component  of  amplitude  A  and  frequency  (Oq.  We  computa¬ 
tionally  achieve  this  variation  via  the  forcing  function 


b  = 


r 


< 


v 


Jt(l+A)  .  f  n  \ 

- - -  sin  I  —  f  if  t  <  ta 

2fa  \ta  J 

-A  (Oo  sin[<Uo  (t  -  tfl)]  if  t  >  ta 


(53) 


We  accelerate  the  flow  from  rest  so  that  U(0)  =0  and  therefore  (53)  yields 


r  i+a 

i  f 71  \ 

1  -  cos  — f 

Va  J . 

u=  < 

.  2 

[  1  +Acos[mo(t-ta)] 

if  t  <  ta , 
if  t  >  ta. 


(54) 


The  first  portion  of  (53)  increases  U  sinusoidally  from  zero  to  some  finite  value,  which  in  this 
study  is  1 +A',  thereupon,  U  is  varied  harmonically.  The  gradual  turning  on  of  the  flow  for 
t  <  ta,  is  necessary  to  ensure  the  continuity  of  the  forcing  function  and  therefore  provide  the 
numerical  stability  of  our  scheme. 

Following  Skopovi  and  Akylas  [23],  we  implement  open  boundary  conditions  (26)  by 
solving  the  governing  equations  (25)  with  the  viscous  terms  added  to  the  right-hand  side  of 
(25a).  The  kinematic  viscosity  (v)  is  spatially  varying  within  the  viscous  layers  surrounding 
the  inviscid  region  of  interest.  Specifically,  near  the  upper  boundary,  we  let 

v  =  0.1M-‘sin(i^f),  (55) 

so  v  increases  gradually  from  zero  at  the  lower  edge  (z  =  zo)  to  v  =  0.1jU_1  at  the  upper 
edge  (z  =  Z)  of  the  viscous  layer.  The  upper  edge  coincides  with  the  top  boundary  of  the 
computational  domain.  In  the  side  layers,  to  ensure  numerical  stability  at  the  inflow  and 
outflow  boundaries,  the  viscosity  is  first  increased  in  a  way  analogous  to  (55),  then  is  held 
constant  at  O.lju-1  for  50  grid  points  and  finally  is  brought  back  to  zero  again  sinusoidally 
over  the  span  of  50  points  (see  Skopovi  [22]  for  details).  In  all  simulations,  we  used  zo=41. 
A  typical  contour  plot  of  the  viscosity  profile  is  presented  in  Fig.  1 1  from  which  it  is  evident 
that  nearly  two  thirds  of  the  computational  domain  are  occupied  by  absorbing  layers.  In  order 
to  accommodate  the  bottom  boundary  of  the  computational  domain  where  we  enforce  the 
inviscid  boundary  condition  (i.e.  zero  velocity  perpendicular  to  the  topography),  viscosity 
in  the  lateral  absorbing  layers  is  linearly  increased  in  the  vertical  direction  from  zero  at  the 
topography  to  the  layer  value  over  20  grid  points,  as  indicated  by  curving  contours  in  the 
bottom  left  and  right  comers  of  Fig.  11. 

Following  Lamb  [16],  we  have  used  a  quadrilateral  grid  with  the  higher  resolution  in  the 
vicinity  of  the  obstacle.  The  grid  is  linearly  stretched  along  both  the  horizontal  and  the  vertical 
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Fig.  11  A  typical  contour  plot  of  the  viscosity  profile.  Viscosity  is  gradually  increased  from  the  inviscid  region  (white)  to  the 
maximum  value  O.ljt-1  (black). 


such  that  the  height  of  the  bottom  cells  is  15  times  smaller  than  the  height  of  the  cells  at  the 
upper  boundary,  while  the  width  of  the  middle  cells  is  15  times  smaller  than  the  width  of 
the  left-  and  right-boundary  cells.  All  simulations  are  carried  out  on  a  numerical  domain  that 
contains  245  cells  in  the  vertical  direction  for  flows  with  Z— 60  and  205  cells  for  flows  with 
Z=50.  In  the  horizontal  direction,  1200  cells  were  used  for  flows  with  Z  =  60  and  860  cells 
for  flows  with  Z= 50  while  the  width  of  the  numerical  domain  is  x  G  [—30, 30]  when  Z  =  60 
and  x  e  [—22, 22]  when  Z = 50.  Lastly,  in  regard  to  temporal  discretization,  upon  impulsively 
accelerating  the  flow  from  rest,  we  use  the  time  step  At  =  0.015. 

The  computational  approach  was  tested  in  three  ways.  We  first  examined  the  performance 
of  viscous  layers  by  comparing  the  behavior  for  the  inviscid  domain  bounded  by  x= ±10  and 
z  =  50  to  that  obtained  with  boundaries  positioned  at  x— ±20  and  z  —  60.  A  comparison  be¬ 
tween  responses  with  £  = 0.6,  jl  =  1/6,  ta  =  0. 05  ji  and  coq  =  12  is  presented  in  Fig.  12.  Results 
suggest  that  even  with  A  =0.75,  layers  at  x  =  ±10  are  effectively  dissipating  disturbances 
by  negligibly  influencing  streamline  patterns  within  the  inviscid  region.  This  comparison  is 
conducted  at  t  =  100,  which  is  the  largest  time  associated  with  our  simulations  involving  the 
varying  free  stream. 

As  the  second  test,  we  confirmed  that  the  computation  in  Fig.  12  reaches  the  steady  state 
and  holds  it  for  a  substantial  amount  of  time.  The  comparison  between  responses  at  r=83  and 
t  =  100  is  shown  in  Fig.  13.  It  is  evident  from  the  figure  that  steady  conditions  are  reached  by 
t  =  83  and  that  the  model  is  holding  this  state  extremely  well. 

As  the  last  check,  we  verified  that  in  the  limit  e  — »0  the  flow  field  with  A =0.1 5  and  a»o = 12 
approaches  the  steady-in-the-mean  solution  predicted  by  the  theory  derived  in  the  previous 
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Fig.  12  Comparison  of  the  response  for  e  =  0.6,  fi  =  1/6,  A  =  0.75,  a»o=  12,  ia  =  0.05/r  and  r  =  100  when  absorbing  layers 
are  positioned  at  x= ±10  (dashed  lines)  with  that  obtained  for  the  same  set  of  parameters  when  the  inviscid  region  extends  to 
jc=±18  (solid  lines). 


Fig.  13  Evolution  of  the  flow  field  for  e=0.6,  fi=  1/6,  A  =0.75,  ta=0.05ji  and  ruo  =12:  dashed  lines /  =  83,  solid  lines  t  =  100. 
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section.  The  agreement  for  e  =  0.2,  fi— 1/6,  and  ta  =  0.05/t  is  presented  in  Fig.  14(a).  With 
an  increase  in  e  to  0.6,  nonlinearity  becomes  important  and  the  two  responses  grow  apart  as 
illustrated  in  14(b). 

We  end  the  discussion  of  unsteady  numerics  by  reporting  that  our  formulation  of  the  dis¬ 
sipative  layer  technique  did  not  work  in  the  absence  of  the  uniform  stream,  that  is,  when 
U  =  cos(fflt)f).  Under  these  circumstances,  our  computations  were  blowing  up.  We  did  not 
explore  the  reasons  behind  this  malfunction  because,  as  pointed  out  by  Lott  and  Teitelbaum 
[18],  the  case  of  reversing  background  flow  is  not  extremely  relevant  to  atmospheric  condi¬ 
tions  that  we  are  concerned  with  in  this  study. 

Gradually  Accelerated  Background  Velocity 

We  explore  the  significance  of  nonlinearity  in  flows  featuring  temporal  variations  in  the 
free  stream  by  gradually  increasing  the  basic  velocity  according  to  (54)  with  A  =  0.  Here,  ta 
then  physically  denotes  the  time  required  to  accelerate  U  from  rest  to  its  terminal  value  U=  1 J 
Due  to  the  fact  that  following  the  initial  acceleration  the  background  flow  is  kept  constant,  we 
expect  the  response  to  approach  the  steady -state  solution  for  U  =  1 ,  which  is  that  of  Long  [  1 7] 
.  With  /l  fixed,  the  problem  under  consideration  is  governed  by  two  parameters,  namely,  £ 
and  ta.  As  it  is  customary  for  a  two-variable  problem,  we  conduct  the  analysis  by  first  holding 
ta  constant  and  varying  £.  Subsequently,  we  fix  £  and  allow  for  changes  in  ta.  In  each  case, 
we  contrast  the  resulting  response  to  that  obtained  with  the  same  value  of  £  and  ta  —  0.05/1, 
which  simulates  the  impulsive  startup. 

Our  fully  nonlinear  simulations  indicate  that  in  the  limit  of  the  vanishing  topography 
height  (f  -» 0),  there  is  a  negligible  difference  between  the  impulsive  and  the  slowly  ac¬ 
celerated  behavior.  This  is  evident  from  the  comparison  of  responses  for  £  =  0.2  and  ju  =  1/6 
presented  in  Fig.  15.  The  contrast  at  t  =  40  in  Fig.  15(a),  an  instant  when  the  gradually  in¬ 
creased  free-stream  velocity  has  just  reached  its  terminal  value,  reveals  negligible  difference 
in  streamline  amplitudes  and  steepness  directly  above  the  barrier.  The  patterns  differ,  how¬ 
ever,  in  the  lee  of  the  obstacle  and  this  discrepancy  is  due  to  transients  of  the  slowly  accel¬ 
erated  response,  which  at  this  time  are  still  present  in  the  domain  of  interest.  They  depart  by 
t  =  60  (Fig.  15(b))  and  after  this  point  the  two  flow  fields  coincide. 

An  increase  in  £  to  0.6  while  keeping  j u  and  ta  constant,  generates  a  profound  difference 
between  the  two  responses.  Particularly,  at  f  =  40,  as  depicted  in  Fig.  16(a),  streamline  am¬ 
plitudes  of  the  slowly  initiated  flow  are  now  nearly  twice  as  large  as  their  rapidly  actuated 
analogues.  Although  in  this  simulation  wave  breaking  did  not  occur,  we  note  that  streamlines 
corresponding  to  the  accelerated  field  are  much  steeper  and  are  nearly  at  the  point  of  over¬ 
turning.  As  transients  move  out,  the  progressively  excited  behavior  once  again  approaches  its 
impulsive  counterpart.  In  comparison  to  the  previously  discussed  case  with  £  =  0.2,  however, 
steady  state  is  reached  at  a  much  later  time,  t  =  140,  as  indicated  in  Fig.  16(b). 

The  described  behavior  for  the  fixed  £  and  two  different  acceleration  times,  0.05jU  and 
40,  suggests  that  for  each  value  of  ta  there  is  a  distinct  nonlinearity  parameter  for  which 
wave  breaking  takes  place.  Correspondingly,  we  have  further  utilized  our  unsteady  numer- 
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times  are  accurate  within  ±2  of  the  plotted  values.  This  tolerance  is  indicated  by  the  error 
bars  surrounding  each  point  on  the  graph.  They  were  obtained  by  refining  the  resolution  of 
the  simulation  for  e= 0.6  until  the  breaking  time  did  not  significantly  change  any  longer.  The 
computed  points  are  also  dependent  on  the  long-wave  parameter  although  for  0<  jj,  <  1  /6  this 
dependency  is  expected  to  be  mild.  This  conclusion  is  based  on  the  fact  that  the  hydrostatic 
solution  wave-breaks  for  £  —  0.85  while  its  counterpart  with  ju  =  1  /6  breaks  at  the  slightly 
lower  value,  £ =0.82. 


Fig.  17  Wave  breaking  as  a  function  of  the  nonlinearity  parameter  and  the  time  to  acceleration  for  fi  =  1/6. 


The  underlying  conclusion  of  this  investigation  is  that  temporal  variations  in  the  incident 
flow  in  combination  with  nonlinearity  may  give  rise  to  large  amplitude  transients,  which  are 
ignored  by  the  previously  constructed  models.  As  a  result,  these  earlier  efforts  underestimate 
the  size  of  the  topography  generated  atmospheric  buoyancy  waves. 

High-Frequency  Sinusoidal  Modulations 

We  now  turn  to  an  investigation  of  the  flow  field  generated  when  the  background  velocity, 
in  accordance  with  (35),  is  a  superposition  of  a  uniform  stream  and  monochromatic  temporal 
variations.  In  discussing  these  flows,  we  limit  our  considerations  to  oscillation  amplitudes  that 
are  smaller  than  the  steady  component  (A  <  1),  so  that  the  basic  velocity  does  not  experience  a 
change  in  direction.  Moreover,  in  this  section,  we  concentrate  on  the  high-frequency  regime 
(coq  -»  oo)  of  these  modulations  and  we  postpone  the  analysis  of  general  frequencies  to  the 
next  section. 
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In  the  limit  Oo  — »  °°,  the  steady -in-the-mean  solution  (47)  for  a  uniform  stratification  and 
H  — >0,  as  analytically  demonstrated  by  Bannon  and  Zehnder  [2],  is  identical  to  the  behavior 
with  the  incident  wind  17  =  1.  In  other  words,  the  transience  in  the  background  velocity  does 
not  impact  the  response;  hence,  at  time  instances  when  the  basic  flow  is  unity,  (47)  coincides 
with  the  linear  solution  of  Long  [17].  Our  numerical  simulations  reveal  that  with  A  =0.75, 
the  two  hydrostatic  streamline  patterns  essentially  coincide  when  (Bo  is  roughly  200.  As  the 
value  of  i u  increases,  however,  the  frequency  at  which  the  agreement  between  the  behavior 
with  A=  0  and  that  with  A  =0.75  takes  place  rapidly  diminishes.  In  particular,  according  to 
our  computations,  for  /i  =  1  /6,  the  accord  is  reached  when  cbq = 6,  which  upon  the  conversion 
to  the  dimensional  form  corresponds  to  the  atmospheric  Brunt-Vaisala  frequency.  We  also 
find  that  aforementioned  conclusions  are  unaffected  by  variations  in  the  buoyancy  frequency 
as  the  steady-in-the  mean  solution  (47)  for  A  =0.75  and  Ct)o  =  6  matches  the  linear  (£  — >0) 
response  of  Davis  [8]  when  the  stratification  is  given  by 

N=  1.5  +  0.5tanh[c(z  —  rf)].  (56) 


with  c  =  2  and  d  =  3. 

The  behavior  exhibited  by  the  linear  response  also  persists  in  the  nonlinear  regime.  This  is 
concretely  evident  from  the  comparison  of  the  resonant  flow  field  in  Fig  2(a)  of  Skopovi  and 
Akylas  [23]  with  that  obtained  for  the  same  stratification  profile  and  the  background  velocity 
(54)  when  ta  =  0.05jl i,  A  =  0.75  and  £Bq  =  12.  The  plot  of  the  two  flow  fields  at  t  =  100, 
illustrated  in  Fig.  18,  reveals  that  streamline  patterns  are  nearly  identical. 

We  lastly  point  out  that  in  agreement  with  conclusions  of  the  previous  section,  we  find  that 
the  amplitude  of  transients  produced  by  temporal  variations  in  the  basic  velocity  increases 
with  the  decrease  in  coq.  This  feature  dominates  at  lower  frequencies  of  temporal  free-stream 
oscillations,  which  we  consider  next. 

General  Frequencies 

In  this  section,  we  discuss  the  flow  field  produced  when  the  modulation  frequency  in 
(35)  is  of  the  order  of  few  hours.  As  before,  we  conduct  the  analysis  by  first  examining 
properties  of  the  linear  steady-in-the-mean  solution  (47)  and  then  we  consider  the  effects  of 
the  nonlinearity  and  transients. 

When  the  frequency  of  monochromatic  modulations  in  (35)  is  coo =0.35  and  jU  =  1/6,  the 
steady-in-the-mean  solution  (47)  qualitatively  retains  the  properties  of  the  response  with  the 
uniform  incident  velocity.  Particularly,  in  the  presence  of  the  resonant  buoyancy  frequency 
given  by 

N2{z )  =  1  +  jU2  {0.25sin(2z)} .  (57) 

the  solution  for  A  =  0.75  exhibits  the  resonant  behavior  in  a  sense  that  the  output  of  our 
numerical  algorithm  does  not  converge  to  a  physically  realistic  solution.  We  were  unable  to 
verify  if  the  resonance  persists  at  lower  values  of  cbo  as  these  simulations  necessitate  a  larger 
range  of  n,  which  in  turn  increases  the  cost  of  our  simulations. 
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Fig.  18  Comparison  between  the  response  with  the  uniform  (dashed  lines)  and  monochromatically  varying  (solid  lines)  basic 
velocity  for  e= 0.6,  fi  =  1/6,  4=0.75,  coo  =  12,  ta— 0.05ft,  F=  1  and  tp=0  at  t  =  100. 

In  addition,  our  computations  indicate  that  the  hydrostatic  limit  of  (47)  for  coq  =  0.35  and 
A  =  0.75,  analogous  to  models  of  Davis  [8]  and  Durran  [9],  can  be  tuned  by  adjusting  the 
tropopause  height.  For  the  stratification  profile  (56)  with  c  =  2  and  <7  =  3,  we  find  that  in  the 
range  d  e  [2,2n]  the  response  is  maximized  when  d  =  3^/2  and  minimized  for  d  =  3jt/4 
(Fig.  19).  These  values  are  in  contrast  to  the  <Mt)  00  limit  where  the  linear  flow  field  is  tuned 
or  detuned  depending  on  whether  d  is  an  even  or  an  odd  multiple  of  jt/2  respectively. 

In  accordance  with  our  findings  earlier,  gradual  monochromatic  temporal  free-stream 
modulations  in  the  presence  of  nonlinearity  produce  transients  whose  amplitude  grows  with 
time.  This  is  concretely  demonstrated  in  Fig.  20,  which  compares  the  response  for  £  =  0.6, 
^t  =  l/6  and  ta  =  0.05/1  at  two  different  times,  namely,  2.5  and  12.5,  when  the  background 
flow  in  Fig.  20(a)  modulates  with  the  amplitude  A  =0.75  and  the  frequency  fflo  =  0.63.  The 
two  times  are  therefore  exactly  one  period  of  the  incident  flow  oscillations  apart.  It  is  evident 
from  Fig.  20(b)  that  streamlines  corresponding  to  t  —  12.5  are  much  larger  in  amplitude  and 
steeper  than  those  associated  with  t  =  2.5.  This  in  turn  suggests  that  similar  to  the  resonant 
buoyancy  frequency  oscillations  studied  earlier,  gradual  variations  in  the  basic  velocity  may 
produce  overturning  of  density  contours  well  below  the  critical  amplitude  of  the  topography 
predicted  by  the  model  of  Long  [17].  Unfortunately,  we  were  unable  to  carry  this  simulation 
to  the  point  of  wave  breaking  as  transients  generated  by  the  gradual  basic  flow  variations 
eventually  become  so  large  that  our  viscous  layers  are  unable  to  absorb  them  even  with  an 
increase  in  the  value  of  v. 
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Fig.  19  Linear  steady-in-the-mean  solution  for  £  —  0.2,  =  1  /6,  A=  0.75  and  =  0.35  when  the  buoyancy  frequency  is  given 

by  (56)  with  c= 2:  d  —  2n/2  (solid  lines)  and  d  =  3lt/4  (dashed  lines). 
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Fig.  20  Transients  produced  by  monochromatic  gradual  temporal  modulations  in  the  free  stream  (a)  background  flow  charac¬ 
terized  by  ta— 0.05/1,  A  =0.75  and  fflo  =  0.63;  (b)  evolution  of  the  response  for  e  =  0.6  and  ji  =  1/6:  dashed  lines  t  —  2.5,  solid 
lines  t=  12.5.  In  this  simulation,  v  =  0.2. 
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